Import

Covid19 Data

Covid19 Japanの公開データ（2020-10-07付）を使います。

　

Area Data

都道府県と地方区分に関するデータは自作データで都道府県コード（pcode）順で各区分を因子化してあります。

　

Data Wrangling

Summarize

オリジナルのデータがどのようになっているかskimrパッケージを用いてサマライズしてみます。元がJSON形式なので、読み込んだ直後は殆どの変量（フィーチャー）が文字型になっているのと欠損値が多いことが分かります。

df %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 88458
Number of columns 23
_______________________
Column type frequency:
character 19
logical 3
numeric 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
patientId 0 1.00 1 8 0 87039 0
dateAnnounced 0 1.00 10 10 0 253 0
gender 9795 0.89 1 1 0 2 0
detectedPrefecture 0 1.00 3 15 0 49 0
patientStatus 84677 0.04 8 23 0 8 0
notes 46068 0.48 1 270 0 39760 1
mhlwPatientNumber 88009 0.01 1 11 0 434 0
prefecturePatientNumber 8324 0.91 5 20 0 80125 0
prefectureSourceURL 57259 0.35 5 224 0 3424 0
residence 17671 0.80 1 38 0 1407 0
sourceURL 586 0.99 1 239 0 7101 0
relatedPatients 79229 0.10 2 259 0 5645 0
knownCluster 86146 0.03 3 88 0 222 0
detectedCityTown 65777 0.26 2 22 0 660 0
cityPrefectureNumber 65842 0.26 1 34 0 22607 2
citySourceURL 76727 0.13 9 317 0 3595 0
deceasedDate 86836 0.02 10 10 0 203 0
deceasedReportedDate 87358 0.01 10 10 0 175 0
deathSourceURL 87452 0.01 14 123 0 618 0

Variable type: logical

skim_variable n_missing complete_rate mean count
confirmedPatient 0 1 0.98 TRU: 87038, FAL: 1420
charterFlightPassenger 88444 0 1.00 TRU: 14
cruisePassengerDisembarked 88447 0 1.00 TRU: 11

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
ageBracket 0 1 33.54 23.17 -1 20 30 50 100 ▃▇▅▂▁

　

Transform

各変量（フィーチャー）を適切な形式に変換し、地域区分でも分析できるように都道府県データと結合します。

x <- df %>% 
  dplyr::mutate(dateAnnounced = lubridate::as_date(dateAnnounced),
                ageBracket = forcats::as_factor(ageBracket),
                gender = forcats::as_factor(gender),
                patientStatus = forcats::as_factor(patientStatus),
                residence = forcats::as_factor(residence),
                cluster = dplyr::if_else(!is.na(knownCluster), TRUE, FALSE)) %>% 
  dplyr::select(dateAnnounced, ageBracket, gender, detectedPrefecture,
                patientStatus, knownCluster, cluster) %>% 
  dplyr::left_join(prefs, by = c("detectedPrefecture" = "pref"))

x
x %>% 
  skimr::skim()
Data summary
Name Piped data
Number of rows 88458
Number of columns 13
_______________________
Column type frequency:
character 2
Date 1
factor 9
logical 1
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
detectedPrefecture 0 1.00 3 15 0 49 0
knownCluster 86146 0.03 3 88 0 222 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
dateAnnounced 0 1 2020-01-15 2020-10-07 2020-08-06 253

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
ageBracket 0 1.00 FALSE 13 20: 21929, 30: 13552, 40: 11018, -1: 9876
gender 9795 0.89 FALSE 2 M: 44322, F: 34341
patientStatus 84677 0.04 FALSE 8 Dec: 1612, Hos: 1266, Hom: 316, Dis: 283
pcode 1029 0.99 FALSE 47 13: 27332, 27: 10978, 14: 7351, 23: 5563
都道府県 1029 0.99 FALSE 47 東京都: 27332, 大阪府: 10978, 神奈川: 7351, 愛知県: 5563
八地方区分 1029 0.99 FALSE 8 関東地: 45816, 近畿地: 17587, 九州地: 10027, 中部地: 9012
広域圏 6124 0.93 FALSE 8 首都圏: 46020, 近畿圏: 17045, 中部圏: 7638, 九州圏: 7292
通俗的区分 1029 0.99 FALSE 11 関東: 45816, 関西: 17045, 東海: 7313, 九州: 7292
fct_pref 1029 0.99 FALSE 47 Tok: 27332, Osa: 10978, Kan: 7351, Aic: 5563

Variable type: logical

skim_variable n_missing complete_rate mean count
cluster 0 1 0.03 FAL: 86146, TRU: 2312

　
文字型を因子型に変換するだけでも、変量（フィーチャー）内の水準の比率が大まかにつかめるようになります。
例えば、年代別の陽性判定者数は20代が最も多く、続いて、30代、40代と働き盛りの世代に多いことが分かります。また、都道府県別では、東京、大阪、神奈川、愛知と成っていますが、地方区分別では、以外にも関東、大阪の次に九州がきていることが分かります。

　

Visualizing

全体傾向

sec_scale <- 50

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n),
                ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)) %>%
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity",
                      fill = "dark green", alpha = 0.5) + 
    ggplot2::geom_line(ggplot2::aes(y = cum / sec_scale),
                       colour = "dark green") +
    ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark red") + 
    ggplot2::scale_y_continuous(
      name = "陽性確定数（日別）・7日間移動平均（濃赤折線）",
      sec.axis = ggplot2::sec_axis(~ . * sec_scale,
                                    name = "累積陽性確定数（折線）")
    ) +
    ggplot2::labs(title = paste0("Covid19, Japan（@", Sys.time(), "）"),
                  subtitle = "", x = "日付")
## Warning: Removed 6 row(s) containing missing values (geom_path).

# plotly::ggplotly()

地方区分別

都道府県別では区分が多すぎるので八地方区分別の積上げ棒グラフを描いてみます。サマリで見たように九州地方が以外にも多いことが読み取れます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
    ggplot2::scale_fill_brewer(palette = "Set1")

折線グラフで表示してみます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) + 
    ggplot2::geom_line(ggplot2::aes(y = n)) +
    ggplot2::scale_colour_brewer(palette = "Set1")

　
さらに分かりやすくするために関東、中部、近畿、九州の四地区を除く他の地区をその他としてまとめてみます。

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
                                                `その他地方` = c("北海道地方",
                                                          "東北地方",
                                                          "中国地方",
                                                          "四国地方"))) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) + 
    ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
    ggplot2::scale_fill_brewer(palette = "Set1")

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
                                                `その他地方` = c("北海道地方",
                                                          "東北地方",
                                                          "中国地方",
                                                          "四国地方"))) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0)) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) + 
    ggplot2::geom_line(ggplot2::aes(y = n)) +
    ggplot2::scale_colour_brewer(palette = "Set1") +
    ggplot2::facet_wrap(~ `八地方区分`, ncol = 2)

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced, `八地方区分`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::pivot_wider(names_from = `八地方区分`, values_from = n) %>% 
  dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>% 
  tidyr::pivot_longer(cols = -dateAnnounced,
                      names_to = "八地方区分", values_to = "n") %>% 
  dplyr::mutate(`八地方区分` = forcats::fct_relevel(`八地方区分`,
                                              "北海道地方", "東北地方",
                                              "関東地方", "中部地方",
                                              "近畿地方", "中国地方",
                                              "四国地方", "九州地方")) %>%
  dplyr::group_by(`八地方区分`) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
                               fill = `八地方区分`)) +
    ggplot2::geom_area(stat = "identity", alpha = 1,
                       position = ggplot2::position_stack(reverse = FALSE)) +
    ggplot2::scale_fill_brewer(palette = "Set3") +
    ggplot2::labs(title = paste0("Covid19（@", Sys.time(), "）"),
                  y = "累計人数")

　
九州が8月頃から急増しているのは、県別に見ると福岡と沖縄での急増が原因と分かります。

x %>% 
  dplyr::filter(`八地方区分` == "九州地方") %>% 
  dplyr::group_by(dateAnnounced, `都道府県`) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::ungroup() %>% 
  tidyr::pivot_wider(names_from = `都道府県`, values_from = n) %>% 
  dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>% 
  tidyr::pivot_longer(cols = -dateAnnounced,
                      names_to = "都道府県", values_to = "n") %>% 
  dplyr::mutate(`都道府県` = forcats::fct_relevel(`都道府県`,
                                              "福岡県", "佐賀県", "長崎県",
                                              "熊本県", "大分県", "宮崎県",
                                              "鹿児島県", "沖縄県")) %>%
  dplyr::group_by(`都道府県`) %>% 
  dplyr::mutate(cum = cumsum(n)) %>% 
  ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
                               fill = `都道府県`)) +
    ggplot2::geom_area(stat = "identity", alpha = 1,
                       position = ggplot2::position_stack(reverse = FALSE)) +
    ggplot2::scale_fill_brewer(palette = "Set3") +
    ggplot2::labs(title = paste0("Covid19, 九州地方（@", Sys.time(), "）"),
                  y = "累計人数")

　

クラスタ比率

x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(`八地方区分`, cluster) %>% 
  dplyr::summarise(n = n()) %>% 
  tidyr::pivot_wider(names_from = cluster, values_from = n) %>% 
  dplyr::mutate(ratio = (`TRUE` / (`TRUE` + `FALSE`) * 100) %>% round(1)) %>% 
  dplyr::rename(`地方` = `八地方区分`,
                `非クラスタ感染者` = `FALSE`, `クラスタ感染者` = `TRUE`,
                `クラスタ比率[%]` = ratio)

Modeling

時系列（TS）分析

日毎の陽性判定者数を時系列データ形式に変換します。判定者数がゼロの日は報告されていないので、最初の陽性判定者が出た日から日日次のカレンダーデータを作成して結合してます。なお、時系列データ形式の周期については日次データなので7日間を周期として設定しておきます（考え方によっては約1ヶ月ごとで4週間=28日を周期にするのもありかと思います）。

ts_x <- x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
  dplyr::select(-date) %>%
  ts(frequency = 1)
  
tsw_x <- x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
  dplyr::select(-date) %>%
  ts(frequency = 7)

tsm_x <- x %>% 
  dplyr::filter(!is.na(fct_pref)) %>% 
  dplyr::group_by(dateAnnounced) %>% 
  dplyr::summarise(n = n()) %>% 
  dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>% 
  dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>% 
  dplyr::select(-date) %>%
  ts(frequency = 28)

時系列データに変換したものをプロットすると可視化の項でプロットした棒グラフと同じ形のグラフになることが分かります。

ts_x %>% 
  plot()

tsw_x %>% 
  plot()

tsm_x %>% 
  plot()

上記からトレンド（長期的傾向）を除いたグラフ。

ts_x %>% 
  base::diff() %>% 
  plot()

tsw_x %>% 
  base::diff() %>% 
  plot()

tsm_x %>% 
  base::diff() %>% 
  plot()

tsw_x %>% 
  stats::decompose() %>% 
  plot()

tsm_x %>% 
  stats::decompose() %>% 
  plot()

スペクトル密度を推定するためのスペクトル密度の推定

ts_x %>% 
  spec.pgram()

tsw_x %>% 
  spec.pgram()

tsm_x %>% 
  spec.pgram()

自己回帰によるスペクトラム解析

spectrum(ts_x, method="ar")

spectrum(tsw_x, method="ar")

spectrum(tsm_x, method="ar")

Infer

ARIMA Model

ARIMA（自己回帰和分移動平均）モデルによる予測を行ってみます。予測に必要なパラメータはステップワイズにより自動的に最適なものが選択されます。

ts_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()

tsw_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()

tsm_x %>% 
  forecast::auto.arima() %>% 
  forecast::forecast() %>% 
  plot()